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Abstract. We investigate a model describing the dynamics of a gas of self-gravitating Brown- 
ian particles. This model can also have applications for the chemotaxis of bacterial populations. 
We focus here on the collapse phase obtained at sufficiently low temperature/energy and on 
the post-collapse regime following the singular time where the central density diverges. Several 
analytical results are illustrated by numerical simulations. 

1. Introduction 

The general study of the static and dynamical properties of a self-gravitating gas 
is a long standing problem in physics. Apart for its clear applications in astrophysics, 
this problem has also many conceptual interests, as for instance, the non-equivalence of 
thermodynamical ensembles for many-body systems with long range interactions |17| . 

In this paper, instead of considering the self-gravitating Newtonian gas (i.e. obeying 
Newton's equations), we study a gas of self-gravitating Brownian particles [TO] subject to 
a friction originating from the presence of an inert gas and to a stochastic force (modeling 
turbulent fluctuations, collisions,...). This system has a rigorous canonical structure where 
the temperature T measures the strength of the stochastic force. Thus, we can precisely 
check the thermodynamical predictions of Kiessling Jl] and Chavanis [S] obtained in the 
canonical ensemble. 

For long-range interacting systems, the mean-field approximation is presumed to be- 
come exact if the thermodynamical limit is taken properly. In a strong friction limit (or 
for large times), the self-gravitating Brownian gas reduces to the Smoluchowski-Poisson 
system, that is a Fokker-Planck equation describing the particle density evolution in the 
presence of a gravitational field $, coupled self-consistently to Poisson's equation, stating 
that the gravitational field is created by the gas itself. 
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This set of equations conserves mass and decreases the Bohzmann free energy |^. They 
describe the competition between the gravitational force which favors a collapsed state 
and the kinetic pressure/diffusion/temperature which tends to spread the particles over 
the entire accessible space. It is thus expected that below a certain critical temperature 
Tc, the system undergoes a situation of "isothermal collapse" which is the canonical 
version of the "gravothermal catastrophe" (TBI . These equations have not been considered 
by astrophysicists because the canonical ensemble is not the correct description of stellar 
systems and usual astrophysical bodies do not experience a friction with a gas (except 
dust particles in the solar nebula Yet, it is clear that the self-gravitating Brownian 
gas model is of considerable conceptual interest in statistical mechanics to understand the 
strange thermodynamics of systems with long-range interactions and the inequivalence 
of statistical ensembles. In addition, it provides one of the first model of stochastic parti- 
cles with long-range interactions, thereby extending the classical Einstein-Smoluchowski 
model to a more general context [Hj. 

In addition, it turns out that the same type of equations occurs in biology in rela- 
tion with the chemotactic aggregation of bacterial populations ^Hl ■ A general model of 
chemotactic aggregation has been proposed by Keller & Segel in the form of two 
coupled partial differential equations. In some approximation, this model reduces to the 
Smoluchowski-Poisson system ^U]. Non-local drift-diffusion equations analogous to the 
Smoluchowski-Poisson system have also been introduced in two-dimensional hydrody- 
namics in relation with the formation of large-scale vortices such as Jupiter's great red 
spot ^1^2 El- These analogies, developed in give further physical interest to our 

Brownian model. 

We now introduce a continuous mass density p{r) in a sphere a radius R, and define 
respectively 

• Total mass: M = J p{r) d^r 

• Energy: E = §MT + i / /9(r)$(r) d^r 

• Entropy: S = f MlnT - / p{r) ln[p(r)] d^r 

At equilibrium, the gravitational potential is given by the Boltzmann-Poisson equation: 



where $ is the gravitational potential, Sd is the surface of the unit D-dimensional sphere, 
and Z is the normalization (partition function). Note that these equations would be the 
same for a Newtonian gas in the mean-field limit, as Eq. ||2Jl simply states that the 
density is given by the Boltzmann weight. The microcanonical version of this model can 
be defined by imposing the total energy E instead of the temperature T. Although the 
equilibrium distributions are the same, the stability limits differ in microcanonical (fixed 
E) and canonical (fixed T) ensembles 5 . A typical phase diagram is displayed in Fig. 
(for D = 3). This illustrates the fact that at low enough temperature in the canonical 
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Figure 1: We show the phase diagram of a self-gravitating Brownian gas in 13 = 3. In 
the canonical ensemble, there is no equilibrium state below the rescaled temperature 
Tc = ryj^ « 0.396, whereas in the microcanonical ensemble, the system collapses below 
the rescaled energy E^. = —Ac « —0.335. The equilibrium phase diagram in all dimensions 
-D, as well as for Langevin particles obeying Tsallis statistics has been studied in [201 17] . 



ensemble, the system does not have any equilibrium state, whereas the same situation 
arises in the microcanonical ensemble at low enough energy. The static phase diagram 
and the detailed stability analysis of the solutions have been extensively studied in PUIIT) 
for isothermal and polytropic distributions and for any dimension of space. 

In this paper, we focus on the dynamical properties of the system when no equilibrium 
state exists. In section 2, after stating the problem, we study scaling collapse solutions 
in the canonical ensemble (for D > 2, including the critical case D = 2) and in the 
microcanonical ensemble (for _D > 2). In both cases, the central density diverges in a 
finite time tcoii (except in D = 2 for T ~ Tc at which tcou — > +00). In section 3, we show 
that the singular state reached at tcoU, which does not coincide with the condensed state 
predicted by statistical mechanics |E1 IHI , is indeed not the final stage of the dynamics. 
We then study the post-collapse stage, which is characterized by the creation of a Dirac 
peak and the existence of a backward scaling dynamics regime. Finally, we illustrate the 
very large time regime where a dilute gas of Brownian particles evolves around a massive 
core. 

2. Collapse dynamics of self-gravitating Brownian particles 
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• Dynamics of the Smoluchowski-Poisson system 

At a given temperature T controlling the diffusion coefficient, the density p{v,t) of a 
system of self-gravitating Brownian particles satisfies the following coupled equations: 

'1 



dt 



^ (VWp + pV$) 



(3) 



A$ SdGp, (4) 

where I? is a diffusion constant. In the present paper, we will consider V — T consistently 
with Einstein's relation. We thus have to solve the Smoluchowski-Poisson system. At 
equilibrium, implying a vanishing current, these equations simply reduce to Eq. H1I2() . 

The more general case where V is a function of p itself can also be of interest (£). For 
example, by taking 'D{p) ~ p^^", one describes a gas of self-gravitating Langevin particles 
displaying anomalous diffusion . The static properties of this system reproduce that of a 
gravitational gas described by Tsallis statistics. The equilibrium distributions correspond 
to polytropes which are associated with the g-entropy Sq — — -^j /(/' — f) dPr dPv , 
where /(r,v) is the phase space density. The parameters q and n are related to each 
other by the relation 

This system can have self-confined equilibrium states depending on n and D. It can also 
undergo gravitational collapse at sufficiently low temperature/energy as illustrated on 
Fig. |5] (lower curve). 

From now on, we take V = T and set M = i? = G = ^ = 1. We shall also restrict 
ourselves to spherically symmetric solutions. The equations of the problem then become 

|^ = v(rvp + pv$), (6) 

A$ = SdP, (7) 

with proper boundary conditions in order to impose a vanishing particle flux on the 
surface of the confining sphere. These read 

^(0,t) = 0, *(1) = ^, t|(1)+p(1) = 0, (8) 

for D > 2. For D = 2, we take $(1) = on the boundary. Integrating Eq. ||7J) once, we 
can rewrite the Smoluchowski-Poisson system in the form of a single integrodifferential 
equation 

9p 1 d ( D-ifrr.dp , p 



r 



The Smoluchowski-Poisson system is also equivalent to a single differential equation 
dM _^/'d^M D-ldM\ 1 ^^dM 
dt \ dr"^ r dr J r^-^ Qr ' 

for the quantity 

M{r,t) - / p{r')SDr'''-Ur\ (11) 
Jo 
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which represents the mass contained within the sphere of radius r. The appropriate 
boundary conditions are 

M(0,t) = iVoW, M(l,t) = l, (12) 

where No{t) = 0, except if the density develops a condensed Dirac peak contribution 
at r = 0, of total mass No{t). It is also convenient to introduce the function s{r,t) = 
M{r,t)/r^ satisfying 

^ -t(— + ^^-^—] + (r— +Ds]s (13) 
dt \9r2 r dr ) \ dr J 

• Self-similar solutions of the Smoluchowski-Poisson system in D > 2 

In ^l|2ni|7], we have shown that in the canonical ensemble (fixed T), the system 
undergoes gravitational collapse below a critical temperature Tc depending on the di- 
mension of space. The density develops a scaling profile, and the central density grows 
and diverges at a finite time tcoii ■ 

We look for self-similar solutions of the form 

Pir,^)=Poit)f{^^), r,-{^Y\ (14) 

where the King's radius tq defines the size of the dense core |2j. In terms of the mass 
profile, we have 

M{r,t) = Mo{t)g(^-^y with Mo{t) ^ p^rl^ , (15) 

and 



g{x) = Sd / f{x')x dx'. (16) 
Jo 

In terms of the function s, we have 

s{r,t)^poit)s(^), with S{x)^^. (17) 



Mt) 

Substituting the ansatz (|17|l into Eq. (|13() . we find that 

^S{x) - P^^xS'{x) ^pli S"{x) + ^^±ls'{x) + xS{x)S'{x) + DS{xf\ , (18) 
dt tq dt \ X J 

where we have set x = r/ro- The variables of position and time separate provided that 
PQ^dpo/dt is a constant that we arbitrarily set equal to 2. After time integration, this 
leads to 

Po{t) = ^{tcoii ~ t)-\ (19) 

so that the central density becomes infinite in a finite time tcoii, which appears as a 
integration constant. The scaling equation now reads 

2S + xS' = S" + ^^^S' + SixS' + DS). (20) 
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Figure 2: In D = 3, we plot s{r,t)/s{0,t) as a function of r/r^it) for different times 
(density range 10^ — 10''). The upper curve corresponds to a constant diffusion coefficient 
V = T, i.e. n = oo, leading to a = 2 jEBEOI- The lower curve corresponds to I? ~ p^/" 
leading to a = ^^3r (the numerical simulations corresponds to n = 4, hence a = 8/3) [7|. 
We compare the numerics to the analytical scaling solutions (dashed lines). 



For D > 2, the scaling solution of Eq. (|20|l was obtained analytically in and reads 

which decays with an exponent a — 2. This leads to (see Fig. |21 upper curve) 

, , 4:{D -2) + D , ^ Ax^ 

= So {D-2 + x^Y^ = D-2 + x^- ^22) 

Note finally that within the core radius tq, the total mass in fact vanishes as t — > tcoii- 
Indeed, from Eq. (|f 5|l . we obtain 

A/(ro(t),i) ~ Po{t)r^{t) ^ T°'\t,oU - t)'"^-\ (23) 

Therefore, the collapse does not create a Dirac peak ("black hole"). 

Near Tc, we find tcou ^ (Tc^T)^^^^ , which is well supported by numerical simulations 
|lUj and by a systematic expansion procedure performed in 0. In addition, the width of 
the scaling regime is St ~ (Tc — Ty^^. This is an estimate of the time tcoU — St from where 
the system enters the scaling regime. Above Tc, the equilibration time characterizing the 
exponential convergence to the stationary solution diverges like r ~ (T — Tc)^^^^ El- 

In [201, we have also studied the collapse dynamics at T = for which we obtained 

Po{t)-^S^\t,ou-t)-\ (24) 
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as previously, but the core radius is not given anymore by the King's radius which vanishes 
for T = 0. Instead, we find 

r-o^Po'/", (25) 

with 



The scaling function S{x) is only known implicitly 

D 

2 



D 



S{x) 



^KxT^S{x), (27) 



where _fir is a known constant (see [3D| for details), 5(0) = -pq:^' ^^^'^ large x asymp- 
totics S{x) ^ f{x) ^ x^". The mass within the core radius is now 

M(ro(t), t) ^ po{t)r^{t) ~ {tcoU - tf^^, (28) 

and it again tends to zero as t — > tcou- Comparing Eq. H23|l and Eq. H28|l suggests that 
if the temperature is very small, an apparent scaling regime corresponding to the T = 
case will hold up to a cross-over time t^, with 

tcoU -U^ (29) 

Above t*, the T ^ scaling ultimately prevails. 

• Scaling and "Pseudo- Scaling" solutions in D ~ 2 

We now consider the critical dimension D = 2 [20]. Above Tc = 1/4, the stationary 
solution can be explicitly computed, and the integrated mass is found to be 

4T 

M(r) = ^— (30) 

Note that M{1) = 1, ensuring that the whole mass is included in this solution. Using 
p = -K^^dM I d{r'^), we find that the density profile is given by 

p{r) ^ — (31) 

with 

ro = V4r- I and p^rl = T. (32) 

Note that the value of and the dependence of and po with the temperature coincide 
with the exact results obtained within conformal field theory 

We now set the temperature to be exactly — 1/4. We then define a{t) in terms of 
the central density 

P(0,t) = ^^. (33) 

TT 

The central density p(0, t) is expected to diverge for t — > +cx), so that a{t) is also expected 
to diverge. Looking for a scaling solution for M{r,t), we find that the scaling function 
coincides with the expression of the stationary solution of Eq. H3U|) and Eq. H31|) . More 
precisely, we find (see Fig. ^ 

^) = 1 f^m 2 + '^{tr'Hr, t), (34) 
1 + a(i)r^ 
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Figure 3: At T = Tc = 1/4, and when the central density has reached the value 
p(0,t) « 1644.8... = {a{t) « 5166.3...), we have plotted the result of the numerical 

calculation compared to our exact scaling form p{r,t) = ^i^i±i(l + a{t)r^)^^ (dashed 
line). The two curves are almost indistinguishable as the relative error is, as predicted, 
of order a^^ ~ 10^''. 



where h{r,t) can also be computed perturbatively in (20). The introduction of the 
next correction to scaling is essential, as it governs the behavior of a{t). Finally, one 
obtains 

S'i^I' + ^C""")]' (35) 

leading to the asymptotic behavior 

a{t) = exp (^^ + V2i^ 1 + 0{t-^/^ \nt) , (36) 

and a similar behavior for the density according to Eq. (|33|l . which is in perfect agreement 
with numerical simulations Note that dXT = Tc the central density does not diverge 
in a finite time tcoii- 

Strictly below Tc, the scaling equation of Eq. p8(l does not have any global physical 
solution. However, we find that by writing 

M{r.t)^- y ^ +M,orir,t), (37) 
ic 1 + a\t)r'^ 

the correction Mcor{r,t) adiabatically satisfies an effective scaling relation of the form 

Meo.(r, t) = a{t)"'^-^hcor{VW», (38) 

where hcor(x) ~ x^^", for large x. Hence, the correction to the density scaling function 
satisfies Pcor{x) ^ x^°'. The index a is a very slowly varying function of time, such 
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that its explicit dependence on time does not alFect the scaling equation for hcor{x) or 
pcor{x) [201 ■ Although the solution is not strictly speaking a true scaling solution, the 
explicit dependence of a on a{t) is so weak that an apparent scaling should be seen 
with an effective a almost constant for a wide range of values of a{t) or density. Hence, 
the total density profile is the sum of the scaling profile obtained at Tc with a T/Tc 
weight (behaving as a Dirac peak of weight T/Tc at t = tcoii, as p{r) ^ r~"'=, with 
ttc = 4 > Z? = 2) and of a pseudo-scaling solution associated to an effective scaling 
exponent slowly converging to a = 2. More explicitly, we find |2U| 



Let us illustrate quantitatively the time dependence of a. For example, Eq. (|39|l respec- 
tively leads to a{a = 10^) = 1.27..., and to a{a = 10^) = 1.34... Thus, for the typical 
values of a{t) accessible numerically of order a ^ 10^, we expect to observe an apparent 
scaling solution with a « 1.3. This is confirmed by the scaling plot of Fig.^J 

Finally, for T < Tc, we note that the central density diverges again in a finite time as 
one has 



• Scaling solutions in the microcanonical ensemble [D > 2) 

A microcanonical ensemble dynamics can be defined formally by considering Eq. 
with a uniform but time dependent temperature T(t) such that the total energy is kept 
strictly constant. The resulting equations increase the Boltzmann entropy and display a 
sort of "gravothcrmal catastrophe" ^01 ■ 

If one looks for a scaling solution, one still has po{t)rQ{t) = T{t), but as the temper- 
ature is not necessarily asymptotically constant near t = tcou , the exponent a character- 
izing the decay of the density scaling function is not determined by simple dimensional 
analysis. In Ref. ^0], we found numerically that the scaling equation Eq. H2U|) has physi- 
cal solutions only for a < amax, with amax — 2.209... for D = 3. We also argued that the 
system will select the exponent amax , since it leads to the maximum increase of entropy. 
This is illustrated on Fig. [3 where the scaling collapse is plotted assuming respectively 
a — 2 and a = Ofmax- Note that if a > 2, the temperature diverges at t = tcou, as 
T{t) ~ po{ty~^^°' . However, in our most recent simulations (and this was confirmed by 
Guerra et al. ^21), we find that as one approaches tcoii, the temperature ultimately sat- 
urates to a finite value implying a = 2. This can be understood by observing that the 
conservation of energy implies that if the temperature diverges, then the potential gravi- 
tational energy should also diverge. However, if a < 5/2 (in = 3), it is straightforward 
to show that the potential energy remains bounded. This exemplifies one of the inherent 
flaws of a model for which the temperature is maintained uniform, as we expect more 




(39) 




(40) 



implying 




(41) 
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Figure 4: At T = Tc/2 = 1/8, we have extracted the next correction to scaling pcor = p — 
T/TcPt=t^- We have then plotted Pco^ir, t)/ Pcor{rma.At)^ t) as a function of a; = r/r^i,^{t), 
where rniax(i) is defined as the location of the maximum of Pcoiif, t). Consistently with the 
apparent scaling observed, we found T'^axlO ^ ^ \/ Pcor(''max(i), i)- For a = 2"^^ x 100 
(n = 1,...,8), we have obtained a convincing data collapse associated to a w 1.3, in 
agreement with the theoretical estimate for a, in this range of a. 



realistically the temperature to grow in the dense core (containing a vanishing mass ac- 
cording to Eq. (|23|l 'l. while remaining finite in the halo. We have shown [7] that in a more 
realistic model where T = T(r, t) is not necessarily uniform, by a phenomenon similar to 
the one leading to solutions for a € [2; amax], the dynamics now selects a unique a > 2. 

As the existence of solutions of Eq. H18(l for a > 2 is interesting and has physical 
implications in a more realistic model with a non uniform temperature, let us mention 
some analytical results about it. In the limit of large dimensions, the scaling equation 
can be perturbatively solved (in powers of D^^). 

Defining z — £^21 (which is of order 0(1)), and 



2 



D + - + 0{D-^), or a;o - Vi^ f 1 + 4^ + 0{D 

zU 



such that >5'(a;o) = ot/ and introducing 

2 D 2{z 



2) 



(42) 



(43) 



' z-1 z{z-l) 

we have obtained the explicit form for the scaling function up to second order in D^^ 



m 



Six) 



a 
D 



a 

Yz 



(44) 
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Figure 5: We plot s{r,t)/s{0,t) as a function of r/ro{t) where ro(t) ^ po{ty^"- We try 
both values a — 2 and a — ctmax = 2.209733..., and compare both data collapses to the 
associated scaling function (dotted lines). The scaling associated to amax is clearly more 
convincing than that for a = 2 (the two sets of curves have been shifted for clarity). 
However, our simulations also suggest that T{t) ^ /3o(^)^~^^""'''' does not diverge at tcoii 
(see the insert where a line of slope 1 — 2/amax ~ 0.09491... has been drawn), so that the 
asymptotic scaling should correspond to a = 2 (see also Guerra et al. [T^l. 



The scaling exponent a is itself a function of S'(O) (or z), defined by 



a 



4 


'1 2 ' 
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15 
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z 
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_z z" z-^ z"' 
This function has a well defined maximum for 
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associated to the value 



amax = 2 + - 



'-Id- 

16 



0{D 



-3^ 



(45) 



(46) 



(47) 



This expansion gives amax — 2.24... in D = 3, in fair agreement with the exact value 
amax = 2.2097... obtained numerically in |10| . In addition, the exponent a = 2 is as- 
sociated to z = DS{0)/2 ~ 2 + 4D^^ + 0{D~^), which coincides with the exact value 
S{0) = 4/(1) - 2) obtained in Eq. PT)). 



3. Post-collapse dynamics in the canonical ensemble 
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• Post-collapse scaling at T = 

According to general results of statistical mechanics E] , the equilibrium state of 
self-gravitating particles in the canonical ensemble is a Dirac peak containing all the mass 
(for D > 2). This is not the structure obtained at tcoii- This implies that the evolution 
must continue in the post-collapse regime. The scenario that we are now exploring (2f| is 
the following. A central Dirac peak containing a mass Nglt) emerges at t > tcoiu whereas 
the density for r > satisfies a scaling relation of the form 

p{r,t)=p^{t)f(^-^, (48) 

where jOo(i) is now decreasing with time (starting from po(i — tcoii) — > -l-oo) and ro(i) 
grows with time (starting from ro(t = tcoii) = 0). As time increases, the residual mass for 
r > is progressively swallowed by the dense core made of particles which have fallen on 
each other. It is the purpose of the rest of this paper to show that this scenario actually 
holds, as well as to obtain analytical and numerical results illustrating this final collapse 
stage. 

For T = 0, the dynamical equation for the integrated mass M{r,t) reads 

dM 1 , 

with boundary conditions 

M(0,t) = iVo(i), M(l,t) = l. (50) 
We define po such that for small r 

M{r,t)~No{t)=po{tf-^ + ... (51) 

Up to the geometrical factor S*^^, po(^) is the central residual density (the residual density 
is defined as the density after the central peak has been subtracted). For r = 0, Eq. (|49l) 
leads to the evolution equation for Nq 

^ = PoiVo. (52) 

As No{t) = ior t < tcoU, and since this equation is a first order differential equation, it 
looks like No{t) should remain zero for t > tcoU as well. However, since poitcou) — +oo, 
there is mathematically speaking no global solution for this equation and non zero values 
for iVo(t) can emerge from Eq. H52(l . as will soon become clear. 
We then define 

4M) (53) 

which satisfies 

ds f ds ^\ Nof ds 
d-t^Vd-r+'^'r^V^Vd-r 



^-+D-^]-'^ + Z^{^-+Ds-po]. (54) 



By definition, we have also s{Q,t) — po{t)/D. 
We now look for a scaling solution of the form 



s{r,t)=po{t)S[-^^], (55) 
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with S'(O) = and 

PoW-roW"", (56) 

where ro is thus defined without ambiguity. Inserting this scahng ansatz in Eq. (|54|l , and 
defining the scahng variable x — r/rg, we find 

J-^ (aS + xS') = S{DS + xS') + -^^{DS + xS' - 1). (57) 
apt) dt por^ x^ 

Imposing scaling, we find that both time dependent coefficients appearing in Eq. (|57f) 
should be in fact constant. We thus define a constant fi such that 

No = fipor^. (58) 

Equation H57|) implies that po ^ (t — tcoii)~^ , which along with Eq. (|56|) implies that 
Nq^ it- Uou)°'°'^'^- We thus find a power law behavior for Nq, which in order to be 
compatible with Eq. (|^ . leads to 



Po(t)= --1 (i-icoH)^'. (59) 



We end up with the scaling equation 

— ^ [aS + xS') + S(DS + xS') + fix-^ {DS + xS' - 1) = 0. (60) 
U — a 

From Eq. H6()|l . we find that the large x asymptotics of S is S{x) ^ a;"". In a short finite 
time after tcou, it is clear that the large distance behavior of the density profile (r 3> tq) 
cannot dramatically change. We deduce that the decay of S should match the behavior 
for time slightly less than tcoU for which S{x) ^ x~^+^. Hence the value of a should 
remain unchanged before and after tcoU ■ Finally, we obtain the following exact behaviors 
for short time after tcoii- 

Po{t) = ^{t-t,Mr\ (61) 



2 

D 



D + 2 

^oW = ( T;V° (t-tcoii)'^, (62) 



No{t) = M (J^^ ' {t-Uou)^. (63) 

We note the remarkable result that the central residual density p(0, t) = S]j^po{t) displays 
a universal behavior just after tcoii, a result already obtained in [2l)|. Moreover, we find 
that No{t) has the same form as the mass found within a sphere of radius ro(t) below 
tcoih given in Eq. 

Moreover, the scaling function S satisfies 

n I 2 / 2D \ 

-S + xS') +S{DS + xS')+px-^{DS + xS' -1)^0. (64) 



\D 

The constant /i is determined by imposing that the large r behavior of s(r, t) (or p(r, t)) 
exactly matches (not simply proportional) that obtained below tcoU, which depends on 
the shape of the initial condition as shown in 20 . Equation (I64|l can be solved implicitly 



14 



C. SIRE AND P.-H. CHAVANIS 



by looking for solutions of the form — z\S{x)]. After cumbersome but straightforward 
calculations, we obtain the implicit form 



1 + —S{x) = 



D 
D + 2 



(65) 



s{x)^^i^ (-2 ) X-—. (66) 



which coincides with the implicit solution given in PO] , Note that S{x) is a function of 
x^ . We check that the above result indeed leads to S{0) — D^^ , and to the large x 
asymptotics 

Note finally that for T = 0, NQ{t) saturates to 1 in a finite time t^nd, corresponding 
to the deterministic collapse of the outer mass shell initially at r = 1. Indeed, using the 
Gauss theorem, the position of a particle initially at r{t = 0) = 1 satisfies 

The position of the outer shell is then 

r{t)^{l-Dty^^, (68) 

which vanishes for tend ~ D~^. 
• Post-collapse scaling at T > 

In the more general case T ^ 0, we will proceed in a very similar way as in the 
previous section. We define again, 



where A^o still satisfies 



We now obtain 



^ = PoNo. (70) 



ds fd^s D + lds\ f ds \ Nq f ds 



T\i^,+—^\ + \^ir+Ds\s + -^{r- + Ds~p,\. (71) 



dt \ dr^ r dr J \ dr J \ d' 

By definition, we have again s{0,t) = po{t)/D. 
We look for a scaling solution of the form 



yr 



s{r,t)^po{t)S[^], (72) 



Mt) 

with S'(O) = D^^. As before, we define the King's radius by 

For t < tcou, we had s{r,t) ~ 4Tr^^ (or S{x) ^ 4a;^^). In a very short time after tcou, 
this property should be preserved, which implies that the post-collapse scaling function 
should also behave as 

Six) - 4a;-2, (74) 
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Figure 6: We plot No{t) for small time (full line). This is compared to 
^^^^^Thcory ^ [l + a{t - tcoii)''] (dashed line), where No{t)'^'^''°'y is given by Eq. ^ with 
fj, = 8.38917147..., and a « 1.7. and b « 0.33 are fitting parameters. Note that the validity 
range of this fit goes well beyond the estimated with — tcoii T^/^ ^ 0.09. The 
bottom insert illustrates the exponential decay of 1 — No{t) ~ e"'*'*. The best fit for A 
leads to A w 5.6362 to be compared to the eigenvalue computed by means of Eq. H83|l . 
A = 5.6361253.... Finally, the top inset illustrates the sensitivity of No{t) to the space 
discretization dx, which introduces an effective cut-off necessary in order to smoothly 
cross the singularity at i = tcoii (a factor 4 in dx between each of the 3 curves). Note the 
small time scale : even the curve corresponding to the coarsest discretization becomes 
indistinguishable from the others for t > 0.448. 



for large x. Inserting the scaling ansatz in Eq. 171|l . we obtain 

(25' + xS') = S" + ^^^S' + S{DS + xS') + -^\{DS + xS' - 1). (75) 
2pq dt X poTq x^ 

Again, this equation should be time independent for scaling to hold, which implies that 
there exists a constant fi such that 

No^fipor^. (76) 

This leads to the universal behavior 



Po{t) 



(f -l)(i-ico»)-^ (77) 
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Figure 7: In the post-collapse regime, we plot p{r,t)/ po{t) as a function of the scaling 
variable x — r/ro{t). A good data collapse is obtained for central residual densities in 
the range 10'^ ~ 10^. This is compared to the numerical scaling functions computed 
from Eq. H78() (dashed line) . The insert shows the comparison between this post-collapse 
scaling function (dashed line) and the scaling function below tcoii which has been rescaled 
to have the same value at a; = 0, preserving the asymptotics: S{x) = (3 + x'^/A)^^ (see 
Eq. |(5lJ; full hue). Note that the post-collapse scaling function is flatter near a; = 0, as 
S{x) - 1/3 ~ x^ (in D = 3) instead of S{x) - 1/3 ~ x^, below tcou- 



We thus end up with the scaling equation 

(2S + xS') + S" + ^^^^S' + S(DS + xS') + nx-^ (DS + xS" - 1) = 0, (78) 
D ~ 2 x 

where fi has to be chosen so that S{x) satisfies the condition of Eq. (|74|l . Its value must be 

determined numerically. Note that for small x, the pre-coUapse scaling function satisfies 

S{x) — S{0) ^ x^, whereas the post-collapse scahng function behaves as 

S{x)-S{0)--x^. (79) 

However, contrary to the T = case, S{x) is not purely a function of x^ . 

Finally, we find that the weight of the central peak has a universal behavior for short 
time after tcoii 

^oit)-f^[-^y^' 'T^/'{t-t..uf^'-\ (80) 
Note that iVo(i) behaves in a very similar manner to the mass within a sphere of radius 
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To below tcoih shown in Eq. 123|) . The behavior of iVo(t) is illustrated in Fig.El while the 
scaling regime is displayed in Fig. [3 

In addition, comparing Eq. H8U|) and Eq. I|63|) . we can define again a post-collapse 
cross-over time between the T 7^ and T — regimes 

U - Uoii ~ T^/', (81) 
which is similar to the definition of Eq. H29|l. 



• Large time limit for T > 

For very large time, that is when almost all the mass has collapsed at r = 0, so that 
No{t) « 1, the residual density satisfies 

p(r,i)~e-^*VW, (82) 

where satisfies the eigenequation 

-A^(r) = T + ^^^') + (83) 

We did not succeed in solving analytically the above eigenequation, and for a given 
temperature, this has to be solved numerically. However, in the limit of very small tem- 
perature, we can apply techniques reminiscent from semiclassical analysis in quantum 
mechanics (T ^ K). We now assume T very small and define h{r) such that 

h{x) dx 



ip{r) = exp 



T 



(84) 



The function h satisfies the following non-linear first order differential equation 

T (^i' + ^^h^ + ^ -h^ = AT, (85) 
with the simple boundary condition 

h{l) = 1. (86) 

In the limit T ^ 0, the term proportional to T in the left-hand side of Eq. H85|l can a 
priori be discarded leading to |21j 

2\Tr^^^ 

Hr) = ====_, (87) 

1 + V 1 - 4ATr2(^-i) 

which is valid for 1 — r ^ T^/"^. Solving perturbatively Eq. (|85() leads to 

^-^ + ^ + - iT-.0), (88) 

where cjj is a D dependent constant. 

In the inverse formal limit of large temperature (although in practice T < Tj,), we 
obtain 

1 

X = D + — -- + ... (T^+00). (89) 

2{D + 2)T ^ ' ^ ' 

The results of Eq. ^ and Eq. ^ are illustrated on Fig. il 



18 



C. SIRE AND P.-H. CHAVANIS 



2.8- 



'h 2.6 
H 2 4 



2.2 



10 



10^ 
10^ 



10 



10 r 



gTTTj 1 I I I llll| 1 I I I llll| 1 I I I llll| 1 I I I llll| 1 I I I 11^ 



10 



JlllJ I I I I 



10' 



10" 




10 10" 
T 



10 



1 



10^ 



_L 



0.1 



0.2 



0.3 



0.4 



Figure 8: We plot A as a function of temperature (insert). The dashed hne is the small 
temperature expansion of Eq. (|88|l . whereas the dotted line is the large temperature 
estimate which is not very accurate in the physically relevant region T < Tc- The main 
plot represents (AT — i) T~^^^ as a function of T^/^ (line and squares), which should 
converge to cd=3 — 2.33810741... according to Eq. JHEl- We find a perfect agreement 
with this value using a quadratic fit (dotted line). 



4. Conclusion 

In this paper, we have illustrated the rich properties of a Brownian model of grav- 
itational collapse in both canonical and microcanonical ensembles. Quantitative results 
have been obtained for any dimension of space D > 2 (including the critical dimension 
D = 2) and any temperature T <Tc (including the peculiar case T = 0). In the micro- 
canonical ensemble, we have shown that although the scaling equation possesses solutions 
corresponding to a > 2, the solution with a = 2 ultimately prevails. However, we have 
argued that in a model for which the temperature is not kept uniform, we should expect 
a value of a greater than 2 to be selected. In the canonical ensemble, we have also shown 
that the singular point t = tcoii is not the final state of the dynamics, which is consis- 
tent with thermodynamical considerations which predicts a totally condensed state at 
r = at equilibrium (14. i5j . We have investigated this post-collapse regime analytically, 
which displays backward scaling solutions. Finally, using semiclassical methods, we have 
described the very large time regime analytically. Note that the post-collapse regime in 
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the microcanonical ensemble is more complex. It is marked by the formation of a central 
body with small mass and small radius but with huge potential energy. This structure is 
reminiscent of a "binary star" in astrophysics. It is surrounded by a very hot halo with 
T +00 that is almost homogeneous. This "binary-halo" structure is the most probable 
structure in the microcanonical ensemble as its entropy S' ^ InT — s- +00 20 . Thus, 
it should be reached in the post collapse regime. However, the Smoluchowski-Poisson 
system becomes ill-defined as T = 00 so that the evolution of the system after tcoii is 
pathological and requires a small-scale regularization [21]. When a small-scale cut-off h 
is introduced, it is found that the mass of the core decreases as h decreases while the 
temperature increases. This corroborates the previous qualitative discussion and gives a 
hint as to how a rigorous description of the post-collapse regime could be undertaken. 

Except for some exact results, most of our analytical results have been obtained by 
perturbative (Eq. (|SH1), Eq. ||Snilv) or non perturbative methods (Eq. (gSIl, Eq. ifTTI) . . . . ) . 
We hope that mathematicians will find the following problems interesting to study in a 
more rigorous way. 

• Concerning the microcanonical collapse scaling equation, it would be interesting to 
justify the existence of the function a[S'(0)] (or S'(0)[q!]), leading to a maximum 
value a = ctmax for physical solutions. 

• The correct mathematical definition of the post-collapse stage following the singular 
point t = tcoii, and the justification of our supposedly exact estimates for po{t) and 
A^o(i) are certainly needed. 
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